# GMM in finance: a tutorial
# authors: Alan De Genaro
# created: 07 / Sept/2020



#####################################################################
# Set working directory 
this.dir <- dirname(parent.frame(2)$ofile)
setwd(this.dir)
######################################################################
# load data
RF.df=read.csv(file="Data_CKLS.csv",sep = ";")

# Load gmm library
library(gmm)
options(digits=4, width=70)

####################################################################
# GMM estimation of just identified ckls model
####################################################################

ckls.moments <- function(parm, x,dt = 1/12)
  {
    # parm = (alpha,beta,sigma,gamma)'
    # data = [r(t+dt)-r(t),r(t)]
    # dt = discretization step
  
  e.hat = as.vector(x[,1] - (parm[1] + parm[2]*x[,2])*dt)
  m2 = e.hat*as.vector(x[,2])
  m3 = e.hat^2 - dt*parm[3]*parm[3]*(as.vector(x[,2])^(2*parm[4]))
  m4 = m3*x[,2]
  g<-cbind(e.hat,m2,m3,m4)
   return(g)
  }
  
data.mat = as.matrix(RF.df)
start.vals = c(0.04,-0.6,sqrt(1.7),1.5)
ckls.mom = ckls.moments(parm=start.vals,x=data.mat)  
names(start.vals) = c("alpha","beta","sigma","gamma")
gmm.ckls.fit = gmm(g=ckls.moments,x=data.mat, start.vals) 

# generates "Tabela 3 (Table 3)" on paper
summary(gmm.ckls.fit)

theta.hat = coef(gmm.ckls.fit)
theta.hat["alpha"]/-theta.hat["beta"]
-theta.hat["beta"]


######################################################################
# Estimate restricted model
######################################################################
cir.moments <- function(parm,x,dt=1/12) {
  # parm = (alpha,beta,sigma)'
  # data = [r(t+dt)-r(t),r(t)]
  # dt = discretization step
  e.hat = as.vector(x[,1] - (parm[1] + parm[2]*x[,2])*dt)
  m2 = e.hat*as.vector(x[,2])
  m3 = e.hat^2 - dt*parm[3]*parm[3]*(as.vector(x[,2]))
  m4 = m3*x[,2]
  g<-cbind(e.hat,m2,m3,m4)
  return(g)
  }

# iterated GMM estimation of overidentified CIR model
start.vals = c(0.06,-0.5,1)
names(start.vals) = c("alpha","beta","sigma")
gmm.cir.fit = gmm(g=cir.moments,x=data.mat, start.vals,type="twoStep",vcov="HAC", optfct="optim")

# generates "Tabela 4 (Table 4)" on paper
summary(gmm.cir.fit)


theta.hat = coef(gmm.cir.fit)
theta.hat["alpha"]/-theta.hat["beta"]
-theta.hat["beta"]
